;+
; NAME:
;       VALIDATE_VELOCITY
;
; PURPOSE:
;       Take the various dense gas tracer spectra (if available), and
;       determine a combined "densegas" spectrum for use with other routines.
;
; CATEGORY:
;       distance-omnibus
;
; CALLING SEQUENCE:
;       VALIDATE_VELOCITY [,csv_file][,VEL_SAVE=vel_save][,HI_SAVE=hisave]
;                         [,/START][,/REAR][,/PLOT][,/VERBOSE] 
;
; INPUTS:
;       NONE
;
; OPTIONAL INPUTS:
;       CSV_FILE  -- Name of the Distance Database CSV file (Default:
;                    'bgps_distance_database.csv') 
;       VEL_SAVE  -- The IDL save file containing the BGPS velocity
;                    structure (Default: ./local/bgps_velocity_struct.sav)
;       HI_SAVE   -- The IDL save file containing the HI sprectrum
;                    structure (Default: ./local/bgps_hi_spectra.sav)
;
;       *Note: (CNUM_LIST) and (START, REAR) are mutually exclusive
;              options: if CNUM_LIST is supplied, then START & REAR
;              are ignored. 
;       CNUM_LIST -- List of BGPS catalog numbers for objects to
;                    validate (as opposed to running the entire
;                    catalog).  (Also, sets /PLOT keyword.)
;       START     -- First BGPS catalog number to process
;                    [Default: #1]
;       REAR      -- Last BGPS catalog number to process [Default:
;                    last entry in the catalog]
;
; KEYWORD PARAMETERS:
;       START   -- First BGPS item to validate  (used for debugging)
;       REAR    -- Last BGPS item to validate  (used for debugging)
;       PLOT    -- Plot to screen (Window 3) the current BGPS item
;       VERBOSE -- Be verbose
;
; OUTPUTS:
;       Adds the combined densegas spectrum to the velocity structure
;       in VEL_SAVE and updates the file on disk.
;
; OPTIONAL OUTPUTS:
;       NONE
;
; COMMON BLOCKS:
;       NONE
;
; EXAMPLE:
;
;
; MODIFICATION HISTORY:
;
;       Created:  09/16/10, TPEB -- Initial Version
;       Modified: 10/08/10, TPEB -- Changing the way the combined
;                                   "densegas" spectrum is obtained,
;                                   add information to Distance
;                                   Database, and code cleanup
;       Modified: 10/14/10, TPEB -- Further refining the combined
;                                   densegas bit -- adding S/N
;                                   spectra, then doing Gaussian fit
;       Modified: 07/14/11, TPEB -- Cleaned up the memory hog by
;                                   undefining the big variables at
;                                   the end of the routine.
;       Modified: 08/26/11, TPEB -- Trying to refine the 'validation'
;                                   of HCO+/N2H+ spectra in terms of
;                                   selecting actual peaks (and not
;                                   noise).
;       Modified: 08/29/11, TPEB -- Added CNUM_LIST optional input for
;                                   only doing validating particular
;                                   objects (used for debugging the
;                                   new 'validation' scheme).
;       Modified: 10/21/11, TPEB -- No longer using interpolated,
;                                   smoothed spectra for velocity
;                                   validation, but rather using the
;                                   original spectra, whose
;                                   filename(s) are now included in
;                                   the velocity structure.
;       Modified: 11/08/11, TPEB -- Removing use of the CSO CS(5-4)
;                                   sprectra from use in the
;                                   repository.
;       Modified: 11/09/11, TPEB -- Updating velocity validation
;                                   scheme to use new routine
;                                   velocity_fit_spectrum.pro to
;                                   return a valid Gaussian fit and
;                                   VLSR based on the original spectra
;                                   (not the smoothed, interpolated
;                                   stuff in the velocity structure).
;       Modified: 10/30/12, TPEB -- Migrated HCO+ / N2H+ spectra over
;                                   to Yancy's (2012) published
;                                   catalog.
;       Modified: 01/31/13, TPEB -- If CS(2-1) is a non-detection,
;                                   then unset the /CS21 has_rv_types
;                                   flag for that object.
;
;-

PRO VALIDATE_VELOCITY, csv_file, VEL_SAVE=vel_save, HI_SAVE=hi_save, $
                       START=start, REAR=rear, CNUM_LIST = cnum_list, $
                       PLOT=plot, VERBOSE=verbose
  
  COMPILE_OPT OBSOLETE
  
  ;; Parse keywords
  IF n_elements(csv_file) EQ 0 THEN csv_file = 'bgps_distance_database.csv'
  IF n_elements(vel_save) EQ 0 THEN vel_save='./local/bgps_velocity_struct.sav'
  IF n_elements(hi_save) EQ 0 THEN hi_save='./local/bgps_hi_spectra.sav'
  plot = KEYWORD_SET( plot )
  verbose = KEYWORD_SET( verbose )
  shift = 0.
  
  nosave = KEYWORD_SET(start) OR KEYWORD_SET(rear) OR KEYWORD_SET(cnum_list) $
           OR KEYWORD_SET(plot)
  
  ;; Get repository params
  defsysv, '!MW', exists = exists
  IF NOT exists THEN galactic_params 
  
  ;; Read in the various pieces of information needed for this routine
  s = read_bgps_csv(csv_file,csv,/ver)
  restore,vel_save,/ver
  restore,hi_save,/ver
  
  ;; If CNUM_LIST supplied, then only do the big loop for those objects
  IF n_elements( cnum_list ) NE 0 THEN BEGIN
     n_bgps   = n_elements( cnum_list ) 
     ind      = WHERE_ARRAY(cnum_list, s.cnum)
     s        = s[ind]
     velocity = velocity[ind]
     start    = 0L
     rear     = n_bgps-1
     plot     = 1b         ;; Since this is for debugging, auto set /PLOT
  ENDIF ELSE BEGIN
     ;; Else, look for START & REAR keywords and adjust accordingly
     n_bgps = csv.nrows
     IF n_elements(start) NE 0 THEN start = long(start) ELSE start = 0L
     IF n_elements(rear) NE 0 THEN rear = long(rear) ELSE rear = n_bgps-1
  ENDELSE  
  
  IF plot THEN window,3, ysize=1000
  
  ;; Loop through each source
  FOR i=start, rear DO BEGIN
     
     ;; Clean up variables used each loop
     undefine,v

     ;; Set variables for this loop
     v          = velocity[i]
     comb       = fltarr(!MW.NVBINS)
     weightnorm = 0.
     this       = s[i]
     CA  = dblarr(3,4)  ;; There are 4 spectral lines
     snr = dblarr(4)
     
     ;; Setup plotting environment  (if desired)
     IF plot THEN BEGIN
        shift = -1.25                          ;; Shift used with /PLOT keyword
        tit = string(v.cnum,s[i].glon_peak,s[i].glat_peak,format=$
                     "('BGPS #',I4,'  l = ',F5.1,'  b = ',F5.2)")
        plot,v.v_std,v.v_std,yr=[-0.25,10.0],/yst,xtit='V!dLSR!n [km s!u-1!n]',$
             ytit='Scaled Intensity',tit=tit,/nodata
     ENDIF
     
     
     ;;======================================================================
     ;; Analyze (& Plot) HCOP
     ;;   We use the full (unwindowed) HCO+ data, which has been
     ;;   baselined by baseline_spec().  To estimate the noise, we use
     ;;   the region outside [-100,200] km/s (no real HCO+ emission).
     ;;   Then, following the prescription in Schlingman et al. (2011)
     ;;   we only use data above 3-sigma to calculate the moments of
     ;;   the distribution -- All using velocity_fit_spectrum.pro.
     
     shift += 1.25
     IF TOTAL(v.hcop) NE 0. THEN BEGIN
        
        yfit = MPFITPEAK(v.v_std, v.hcop, A, NTERMS=3, /POSIT)
        CA[*,0] = A
        
        comb       += v.hcop / (v.sig_hcop * v.sig_hcop)
        weightnorm += 1.    / (v.sig_hcop * v.sig_hcop)
        
        IF verbose THEN print,'HCOP:  ',A
        IF plot THEN al_legend,position=[190,shift+0.4],['**'],box=0
     ENDIF ELSE yfit = fltarr(!MW.NVBINS)
     
     IF plot THEN BEGIN
        oplot,v.v_std,v.hcop   / 1. + shift, thick=1
        oplot,v.v_std,yfit / 1. + shift, color='c0ff00'x, thick=2
        al_legend,position=[-90,shift + 0.4],['HCO!u+!n'],box=0
     ENDIF
     
     
     ;;======================================================================
     ;; Analyze & Plot N2HP
     ;;   We use the full (unwindowed) N2H+ data, which has been
     ;;   baselined by baseline_spec().  To estimate the noise, we use
     ;;   the region outside [-100,200] km/s (no real N2H+ emission).
     ;;   Then, following the prescription in Schlingman et al. (2011)
     ;;   we only use data above 3-sigma to calculate the moments of
     ;;   the distribution -- All using velocity_fit_spectrum.pro.
          
     shift += 1.25
     IF TOTAL(v.n2hp) NE 0. THEN BEGIN
        
        yfit = MPFITPEAK(v.v_std, v.n2hp, A, NTERMS=3, /POSIT)
        CA[*,1] = A
        
        comb       += v.n2hp / (v.sig_n2hp * v.sig_n2hp)
        weightnorm += 1.    / (v.sig_n2hp * v.sig_n2hp)
        
        IF verbose THEN print,'N2HP:  ',A
        IF plot THEN al_legend,position=[190,shift+0.4],['**'],box=0
     ENDIF ELSE yfit = fltarr(!MW.NVBINS)
     
     IF plot THEN BEGIN
        oplot,v.v_std,v.n2hp   / 1. + shift, thick=1
        oplot,v.v_std,yfit / 1. + shift, color='c0ff00'x, thick=2
        al_legend,position=[-90,shift + 0.4],['N!d2!nH!u+!n'],box=0
     ENDIF
      
     
     ;;======================================================================
     ;; Analyze & Plot CS21
     ;;   CS21 data from ARO 12m.  
     
     shift += 1.25

     yfit = velocity_fit_spectrum(v, VLSR=vlsr_cs21, /CS21, IS_SPEC=is_spec, $
                                  FA=A, SN=sn, VERBOSE=verbose)
     CA[*,2] = A
     snr[2]  = sn
     IF is_spec THEN BEGIN
        IF plot AND sn GT 0 THEN $
           al_legend,position=[190,shift+0.4],['**'],box=0
        
        comb       += yfit / (v.sig_cs21 * v.sig_cs21)
        weightnorm += 1.   / (v.sig_cs21 * v.sig_cs21)
        sig = v.sig_cs21
     ENDIF ELSE BEGIN
        sig = max(v.cs21)
        v.has_rv_types = v.has_rv_types AND (2^12 - 5)
     ENDELSE
     
     ;; Add information to the Distance Database (not yet implimented...)
     
     ;; Plot if desired
     IF plot THEN BEGIN
        oplot,v.v_std,v.cs21  / 2. + shift, thick=1
        oplot,v.v_std,yfit / 2. + shift, color='c0ff00'x, thick=2
        al_legend,position=[-90,shift + 0.4],['CS'],box=0
        al_legend,position=[-70,shift + 0.4],$
                  ['SN = '+string(sn,format="(F0.1)")],box=0
     ENDIF
     
     
     ;;======================================================================
     ;; Analyze & Plot NH3
     ;;   NH3 data from GBT.  This data a is little funny, since the
     ;;   inversion transitions of ammonia have fine structure lines
     ;;   surrounding the main peak.  However, I have the fits from
     ;;   Miranda Dunham (using the full QM fits to the HFS) for
     ;;   v_lsr, tmb, and linewidth.  So, in the
     ;;   'import_densegas_spectra.pro' routine, I create these
     ;;   Gaussians.  So, here, it's easy!
     
     shift += 1.25
     IF TOTAL(v.nh3) NE 0. THEN BEGIN
        
        yfit = MPFITPEAK(v.v_std, v.nh3, A, NTERMS=3, /POSIT)
        CA[*,3] = A
        
        comb       += v.nh3 / (v.sig_nh3 * v.sig_nh3)
        weightnorm += 1.    / (v.sig_nh3 * v.sig_nh3)
        
        IF verbose THEN print,'NH3:  ',A
        IF plot THEN al_legend,position=[190,shift+0.4],['**'],box=0
     ENDIF ELSE yfit = fltarr(!MW.NVBINS)
     
     IF plot THEN BEGIN
        oplot,v.v_std,v.nh3   / 1. + shift, thick=1
        oplot,v.v_std,yfit / 1. + shift, color='c0ff00'x, thick=2
        al_legend,position=[-90,shift + 0.4],['NH!d3!n'],box=0
     ENDIF
     
     
     ;;======================================================================
     ;; Plot combined spectra
     ;;   Plot the combined spectrum made from the product of dense
     ;;   gas spectra that actually exist.
     
     shift += 1.25
     
     IF TOTAL(CA[0,*]) NE 0. THEN BEGIN
        
        ;; To create the combined dense gas spectrum, we use the
        ;; Gaussian fit parameters
        MA = CA[*,WHERE(CA[0,*] NE 0., nsp)]
        
        ;; Check for "valid" spectra that have offset VLSR > 8 km/s
        IF (max(MA[1,*]) - min(MA[1,*])) GE 8. THEN BEGIN
           
           ;; If NH3 present, use it exclusively and move on
           IF CA[0,3] NE 0. THEN FA = CA[*,3] ELSE BEGIN
              ;; Otherwise, use line with largest SNR
              maxsn = max(snr,maxi)
              FA = CA[*,maxi]
           ENDELSE
           
        ENDIF ELSE BEGIN
           ;; Then, work with everyone who DID line up nicely
           IF nsp GT 1 THEN FA = median(MA, DIMENSION=2) ELSE $
              FA = MA
        ENDELSE
        
        ;; Set peak = 1., since arbitrary, anyways.
        FA[0] = 1.d
                
        yfit = gauss_1(v.v_std, FA)
        v_densegas = FA[1]
        
        IF verbose THEN BEGIN
           print,string(fA,format="('DenG:  ',F5.2,'  ',F6.2,'  ',F5.2)")
           print,' '
        ENDIF
        IF plot THEN al_legend,position=[190,shift+0.4],['**'],box=0
     ENDIF ELSE BEGIN
        max=1.
        v_densegas = -1000.
        yfit = fltarr(!MW.NVBINS)
     ENDELSE
     
     IF plot THEN BEGIN
        oplot,v.v_std,comb  /max + shift, thick=1
        oplot,v.v_std,yfit /max + shift, color='c0ff00'x, thick=2
        oplot,v.v_std,comb*0.+0.3 + shift,color='ab60f6'x,linestyle=2
        al_legend,position=[-90,shift + 0.4],['Dense Gas'],box=0
        ;;print,'DG_FIT: ',FA
     ENDIF 
     
     ;; Store the Gaussian Fit and combined spectrum in the velocity structure
     velocity[i].densegas = yfit
     velocity[i].vlsr     = v_densegas
     
     
     ;;======================================================================
     ;; Analyze & Plot GRS 13CO
     ;;   13CO sometimes has multiple peaks (even with "on-off"
     ;;   technique), and has small fractional noise, so plot
     ;;   directly, normalized to peak "on-off" intensity.
     
     shift += 1.25
     IF TOTAL(v.grs13co) NE 0 THEN BEGIN
        max = max(v.grs13co)
        ;; comb *= v.grs13co
        IF plot THEN al_legend,position=[190,shift+0.4],['**'],box=0
     ENDIF ELSE max = 1
     
     IF plot THEN BEGIN
        oplot,v.v_std,v.grs13co/max + shift
        al_legend,position=[-90,shift + 0.4],['GRS !u13!nCO'],box=0
     ENDIF
     
     
     ;;======================================================================
     ;; Plot HISA spectrum
     ;;   For comparison, plot the HISA spectrum above the combined
     ;;   spectrum to see how everything fits together.
     
     shift += 1.25
     IF TOTAL(hi[i].spectrum) NE 0 THEN BEGIN
        max = max(hi[i].spectrum)
        IF plot THEN al_legend,position=[190,shift+0.4],['**'],box=0
     ENDIF ELSE max = 1
     
     IF plot THEN BEGIN
        oplot,hi[i].v_std,hi[i].spectrum/max + shift,color='0000ff'x
        al_legend,position=[-90,shift + 0.4],['HISA'],box=0
     ENDIF
     
     
     ;;======================================================================
     ;; If such an entry exists in the Distance Database, overplot as
     ;; a vertical line the HHT velocity as quoted in the July24.csv
     ;; file from Wayne.
     
     ;; IF s[i].rv1_type NE 0 AND plot THEN vline,s[i].rv1_mag,color='00ff00'x
     IF v_densegas NE -1000 AND plot THEN vline,v_densegas,color='00ffff'x
     
     ;;======================================================================
     ;; Plot the Tangent Velocity for comparison
     
     IF plot THEN vline,!MW.V0*(1. - sin(s[i].glon_peak * !dtor))
     
     
     IF plot THEN $
        IF n_elements(cnum_list) NE 0 THEN res = GET_KBRD() ;ELSE wait,1.0
     
  ENDFOR
  
  ;; Check to only save the velocity structure if doing the whole
  ;; thing (no plotting, start, rear, or cnum_list)
  IF ~ nosave THEN $
     save,velocity,filename='./local/bgps_velocity_struct.sav',/ver
  
  ;; Clean up the memory, dude
  undefine,velocity
  undefine,grs
  undefine,hi
  
END
